Vortices in layered superconductors with columnar pins: A density functional study 
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We use numerical minimization of a model free energy functional to study the effects of columnar 
pinning centers on the structure and thermodynamics of a system of pancake vortices in the mixed 
phase of highly anisotropic layered superconductors. The magnetic field and the columnar pins are 
assumed to be perpendicular to the layers. Our methods allow us to study in detail the density 
distribution of vortices in real space. We present results for the dependence of the average number 
of vortices trapped at a pinning center on temperature and pinning strength, and for the effective 
interaction between nearby pinned vortices arising from short-range correlations in the vortex liquid. 
For a commensurate, periodic array of pinning centers, we find a line of first order vortex lattice 
melting transitions in the temperature T vs. pin concentration c plane, which terminates at an 
experimentally accessible critical point as c is increased. Beyond this point, the transition is replaced 
by a crossover. Our results should also apply, with little change, to thin-film superconductors with 
strong point pinning. 



I. INTRODUCTION 

In the mixed phase of type-II superconductors, mag- 
netic flux penetrates the sample as quantized vortex lines 
which form a special physical system known as "vor- 
tex matter" . The fascinating equilibrium and dynamical 
properties of vortex matter in the mixed phase of high- 
temperature superconductors (HTSCs) have prompted 
considerable experimental and theoretical attention^ for 
more than a decade. Because of enhanced thermal fluc- 
tuations, the Abrikosov lattice in very pure samples of 
these highly anisotropic, layered materials is observed to 
undergo a first order melting transitionEI into a resistive 
vortex liquid as the temperature is increased. 

Thepproperties of the mixed phase of HTSCs are also 
knownS to be strongly affected by the presence of pinning 
centers, either intrinsic to the sample or artificially gen- 
erated. Understanding the effects of pinning in these sys- 
tems is very important for practical applications because 
the presence of pinning strongly influences the value of 
the critical current in the mixed phase. Columnar pin- 
ning arising from damage tracks produced by heavy-ion 
bombardment has received much attention in this con- 
text because such extended defects parallel to the direc- 
tion of the average magnetic flux are highly effective&tj in 
increasing the critical current by localizing vortex lines 
along their length. Columnar defects produce "strong 
pinning" in the sense that the pinning potential of a 
defect is sufficiently strong to pin a vortex line at low 
temperatures. Heavy-ion irradiation generally produces 
a random array of parallel columnar defects. The ef- 
fects of such an array of extended defects on the proper- 
ties of the mixed phase of HTSCs have been extensively 



studied-.experimentallyB~ theoreticallyaO and numeri- 
callyBEJ. The thermodynamics of a collection of vortex 
lines in the presence of a parallel array of random colum- 
nar pins has been analyzecfla by mapping the problem to 
the quantum mechanics of a system of two-dimensional 
interacting bosons in an external random potential. The 
main prediction of such theoriesjSpthe existence of a low- 
temperature "Bose glass" phaseotj, separated by a con- 
tinuous phase transition from a high-temperature, entan- 
gled liquid of vortex lines. The theoretically predicted 
scaling hehavior near the Bose glass transition has been 
observed!! A random array of columnar pins also affects 
the equilibrium properties of the high-temperaturei-wor- 
tex liquid, leading to the occurrence of anomalicsliil in 
the reversible magnetization curve near B = B^, where 
Bcj, = Pp^a {p P is the areal density of columnar pins and 
$0 = hc/2e is the flux quantum) is the so-called "match- 
ing field", and B is the magnetic induction that deter- 
mines the areal density po of vortex lines (po = B/Qq). 

It is also possible, through .the use of a variety of 
nanofabrication techniqucsEf1I3, to create periodic ar- 
rays of strong (in the above mentioned sense) pinning 
centers in thin-film superconductors. The interplay be- 
tween the lattice constant of the pin array (determined 
by B^) and the intervortex separation (determined by B) 
can lead to a variety of interesting commensurability ef- 
fects in such systems. Some of these effects have been_ob=. 
served in recent experiments. Imaging experimentsEj'EZl 
have shown the formation of ordered structures of the 
vortex system at low temperatures for commeasurate val- 
ues of B/B^. Magnetization measurementstHl in the ir- 
reversible (vortex solid) regime have demonstrated the 
occurrence of anomalies at harmonics of B^. The effec- 
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tiveness of pinning at integral values of B / B^ has been 
foundtfl to produce regularly spaced sharp minima in the 
resistivity versus field curve. A pinning-induced recon- 
figuration o£,the vortex lattice has been observed in an 
experimentll3 on a thin-film superconductor with a rect- 
angular array of pinning centers. Some of these effects 
have been studied theoretically, using analytidl3 and nu- 
mericaO methods. Bulk HTSC samples with periodic 
arrays of columnar pins have not been fabricated yet, 
but the technology for doing this appears to be within 
reachEJ. 

A periodic array of strong pinning centers should have 
significant effects on the melting transition of the vortex 
lattice. We consider here the situation where B > B^, 
that is, when the pin array is relatively dilute. If, in addi- 
tion, the value of B is such that the melting temperature 
of the vortex lattice in the pure system is substantially 
lower than the superconducting transition temperature in 
zero field, then each pinning center would pin a vortex at 
temperatures comparable to the melting temperature of 
the pure vortex lattice. However, the interstitial vortices, 
which would be present whenever the number of pinning 
centers is smaller than the number of vortices (assuming 
that each defect can pin at most one vortex), may un- 
dergo a sharp melting transition. This would certainly be 
the case in the limit where the spacing of the pin array is 
sufficiently large. Since the vortices pinned at the defects 
produce a periodic potential for the interstitial ones, the 
melting transition of the latter is an example of a solid to 
liquid transition in the presence of an external periodic 
potential. Evidence for such melting-xif interstitial vor- 
tices has been found in experimentsE3liia on thin-film su- 
perconductors with periodic pinning. However, the ther- 
modynamic behavior at the melting transition has not 
been characterized in these experiments. The effects of 
a periodic potential on the melting of two-dimensional 
solids have-been investigated earlier using analytictll and 
numerical23'E2l methods. We are not aware of any theo- 
retical study of the effects of a periodic array of columnar 
pins on the vortex lattice melting transition in three di- 
mensions. 

In this paper, we report the results of a study of 
the equilibrium properties of the mixed phase of highly 
anisotropic, layered superconductors in the presence of 
columnar pins. We consider a geometry in which both 
the magnetic field and the direction of the columnar 
defects are perpendicular to the superconducting lay- 
ers, fim^-study is based on a model free energy func- 



tionaEJTfd for a system of "pancake" vortices lying on 
the superconducting layers. We consider the limiting 
case of a vanishingly small Joscphson coupling between 
the layers, so that the pancake vortices on different lay- 
ers interact via only theit electromagnetic coupling. As 
shown in earlier studiesE3~E3, this limit is appropriate 
for describing extremely anisotropic Bi- and Tl-based 
HTSCs. Xhe Ramakrishnan-Yussouff (RY) free-energy 
functional^.! used in the-ipesent work is the same as that 
used in earlier studiesc3i2a of vortex lattice melting in 



pure systems. The same free energy functional was also 
used, in combination with the,ceplica method for treating 
quenched disorder, in a studjC3 of the effects of random 
point pinning on the melting line in the B — T plane. 
In these earlier studies, the density distribution in the 
crystalline state was expressed in terms of a few "order 
parameters" and the free energy was minimized with re- 
spect to these parameters. In the present work, we use 
a different method which is more powerful and more ap- 
propriate for describing in detail pinning-induced inho- 
mogeneities of-the local density. This method, developed 
in our studiesE3 of the hard-sphere system, involves di- 
rect numerical minimization of a discretized version of 
the free energy functional. Since both the magnetic field 
and the direction of the columnar pins are assumed to 
be perpendicular to the layers, the time- averaged local 
density of pancake vortices must be the same on all the 
layers. This simplification makes the problem effectively 
two-dimensional and allows a high-precision numerical 
investigation of the effects of columnar pins on the struc- 
ture and thermodynamics of the vortex system. Further- 
more, our results should apply, with little change, to thin 
film superconductors with strong pinning. 

The model considered in our work is defined in Sec- 
tion ||, where the method of calculation is also described. 
We then consider (subsection III A ) the effects of an iso- 
lated columnar pin on the structure of the vortex liquid 
in the vicinity of the pin. This is done mainly for test- 
ing the systematics of our numerical method and also 
for determining appropriate values of the pinning poten- 
tial to be used in subsequent calculations. We choose 
throughout the computations numerical parameter val- 
ues appropriate for BSCCO. We determine the suitable 
choice of the discretization scale in order that our numer- 
ical method provides an accurate account of the density 
inhomogeneities produced by the trapping of a vortex at 
a pinning center. We also determine the range of val- 
ues of the pinning potential strength for which nearly 
one vortex is trapped at a pinning center in the tem- 
perature range of interest. The strength of the pinning 
potential is kept fixed in this range in our subsequent 
work: pinning of multiple vortices at a pinning center is 
not considered because this is rarely observed in experi- 
ments. Next, in subsection IIIB, we consider the effects 



of two neighboring pinning centers on the liquid-state 
properties. An "effective interaction" between vortices 
trapped at the two pinning centers is obtained by cal- 
culating the free energy as a function of the separation 
between the pinning centers. This effective interaction 
is found to oscillate with distance. This study and the 
one-pin calculation mentioned above complement, in a 
sense, the analytic work of Ref. |l8| where the RY free- 
energy functional was used to analyze the structure and 
magnetization of a two-dimensional vortex liquid in the 
presence of strong pinning. However, we consider here a 
three-dimensional system with columnar pins, instead of 
a two-dimensional system with point pinning as in Ref. 
[Ts| . Also, the numerical direct minimization method used 
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in the present work is more accurate than the analytic 
variational method in the earli er stu dy. 

We next study (subsection IIIC) the freezing of the 
vortex liquid in the pure system. This is done primar- 
ily for checkiug.the method against the results of earlier 
calculationscJcJ. We find results in excellent agreement 
with those of earlier studies. Our new method also pro- 
vides a very detailed and accurate account of the distri- 
bution of the densi ty nea r a lattice point. We then pro- 
ceed, in subsection HID, to consider the melting transi- 
tion of interstitial vortices in a commensurate, triangular 
array of columnar pins. As discussed above, this transi- 
tion provides a physical realization of three-dimensional 
melting in the presence of a commensurate periodic po- 
tential. Defining the concentration c of pinning centers as 
c = B^/B, we consider values of c given by l/l 2 where I 
is an integer. For small concentrations of pinning centers 
(I > 6), we find a first-order melting transition from a 
crystalline solid to an inhomogeneous liquid. As the pin 
concentration is increased, the transition temperature in- 
creases and the latent heat and the jump in the crys- 
talline order parameter at the transition decrease. We 
find that this line of first-order transitions terminates at 
a critical point beyond which the thermodynamic transi- 
tion is replaced by a sharp crossover. This critical point 
is a rare, experimentally realizable example of continuous 
melting in three dimensions. We show that a simple Lan- 
dau theory provides a semi-quantitative understanding of 
most of our results. Some of our most salient results on 
the melting transition in the presence of a—periodic pin 
array were summarized in a recent Lettered. Here, we 
present many details which could not be included in that 
short paper. Section [V contains a summary of the main 



results and some concluding remarks. 



II. MODEL AND METHODS 

As explained in the introduction, we perform in this 
paper a numerical study using density functional theory, 
which involves, as its foundation , a model free energy 
functional appropriate to a system of pancake vortices in 
a layered superconductor . Density functional methods 
have long been usedE-El with great success in the study 
of solidification phenomena in ordinary fluid systems. Al- 
though the theory is basically mean-field based, it works 
very accurately in the description of first order melting. 
We have ourselves performed extensive numerical studies 
of a density vs. disorder strength phase diagram of aJaard 
sphere system in the presence of quenched disordered us- 
ing a methodology quite similar to that employed in the 
present work. D en s ity functional methods have also been 
successfully usedt3t3 to study the melting of the vortex 
lattice in layered superconductors without pinning. 

The starting point of our calculation is an expression 
for the free energy of the system, written as a functional 
of the time averaged local density. In our case the rel- 



evant density is /?(i,r), the time averaged areal density 
of pancake vortices at point r on the ith layer. In the 
homogeneous vortex liquid state in the absence of pin- 
ning, this density is uniform and it is given in terms of 
the magnetic induction B by po = B/&o where <3?o is 
the superconducting flux quantum. It is customary and 
convenient to introduce a length oq through the relation 
ttciqPo — L We will use oo as our standard unit of length 
in terms of which other lengths will be given and we will 
usually also normalize densities in terms of po- 
We write the free energy functional in the form: 



F[p}=F RY [p}+F p [p}. 



(2.1) 



The first term in the right-hand side of Eq.(2.1) is the free 
energy of the vortex system in the absence of pinning, 
while the second includes the pinning effects. Since the 
potential produced by a collection of straight columnar 
pins perpendicular to the layers is the same on every 
layer, the time-averaged density of vortices at any point 
r must be the same on all layers, i.e. p(i, r) = p(r) for all 
i. Then, the free energy per layer corresponding to the. 
first term in the right side of Eq.(2.1) may be writtcnE.il 
in an effectively two-dimensional form: 



pF RY \p]= / dr{p(r)ln(p(r)/p )-6p(v)} 



(1/2) / dr / dr'CQr - r'\)5p(r)Sp{r'), 



(2.2) 



where f3 is the inverse temperature. We have defined 
Sp(r) = p(r) — p as the deviation of p(r) from p , the 
density of the uniform liquid and taken our zero of the 
free energy at its uniform liquid value. The function C{r) 
is a static correlation function that contains all the re- 
quired information about the interactions in the system. 
It is given by Y] C(n, r) where n is the label denoting 
layer separation, r is the the in-plane .distance and C(n, r) 
is the direct pair correlation functional of a layered liquid 
of pancake vortices with areal density po. 

Strongly anisotropic layered superconductors can be 
described in terms of the Lawrence-DoniachL 2 ] Hamilto- 
nian. Using that Hamiltonuuji as the starting point, it 
is relatively straightforwardcJ to calculate C(n,r) and 
C(rl-|Using the hypernetted chain (H NC) approxima- 
tionEj. Since the interaction term in Eq. (pj) is of a con- 
volution form, it is numerically most efficient to deal with 
it in wavevector space. This and the use of fast Fourier 
transform (FFT) methods reduces the computation of 
the interaction term in the free energy to a single sum. 
One begins with the expression for the Fourier transform 
of v(n, r), the two-body vortex-vortex interaction, which, 
in the limit of vanishingly .small Josephson coupling be- 
tween the layers, is givenocj by: 



Mk) 



2nTX 2 [kl + (4/d 2 ) sin 2 (fc z d/2)] 
k 2 ± [l + A 2 fc2_ + 4(A 2 /d 2 ) sin 2 (fc z d/2)] ' 



(2.3) 



where k z and k± are respectively the components of k 
perpendicular and parallel to the layer plane, d is the 
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layer spacing, and A(T) the penetration depth in the layer 
plane. The dimensionless quantity T which determines 
the strength of the interactions is givenE3 by: 

r = /3d^/87r 2 A 2 . (2.4) 

In coordinate space, u(0, r) is repulsive and logarithmic 
in r while v(n ^ 0,r) is also logarithmic, but attractive 
and weaker than the intralayer potential by a factor of 
roughly (d/\)e~ nd / x . Beginning with an interaction of 
this form, the HNC procedure of Ref. ^6] can be used to 
numerically compute C(k) for the appropriate values of 
the relevant parameters. The quantity C(k±), the two- 
dimensional Fourier transform of C (r) , is then obtained 
by setting k z — in C(k). 

The second term in Eq.(|]l]) represents the contribu- 
tion of pinning to the free energy per layer. It is given 
by: 

l3F p [p] = J drV p (r)6p(r) (2.5) 

where V p (r) is the dimensionless (normalized by fcgT) 
pinning potential at point r. This quantity can be written 
as V p (r) — J2j v p(\ r ~ Rj | ) j where the sum is over all 
pinning centers located at the points {Rj} on a plane, 
and v p (r) is the dimensionless form of the potential at 
r due to a pinning center at the origin. We rt4ke this 
potential to be of the truncated parabolic forrr£3: 

Vpir) = -V [l - (r/r Q ) 2 ]0(r - r), (2.6) 

where ro is the range of the pinning potential. We will 
write the dimensionless strength parameter Vq as Vq — 
aT and the quantity a will be chosen, as explained below, 
so that the pinning is strong enough to localize one vortex 
at a pinning center at the temperatures of interest, but 
not so strong that more than one vortex is bound to a 
pinning center. 

In order to carry out numerical work, we have to dis- 
cretize our system. We introduce for this purpose a com- 
putational triangular lattice of size L. On the sites of 
this lattice we define density variables pi = p(vi)v, where 
p(vi) is the density at mesh point i and v the area of 
the unit cell in the computational lattice, proportional 
to the square of its lattice constant h. We have L = Nh, 
so that the computational lattice has iV 2 sites. Periodic 
boundary conditions are used in all our calculations. 

Our basic procedure is to minimize the free energy of 
the system given certain values of the relevant parame- 
ters and the appropriate initial conditions, that is, some 
initial set of values for the computational variable pi. 
Finding the minima of the free energy is not at all triv- 
ial, since one is minimizing a function of a very large 
number of variables (we have used values of N up to 
2048 as discussed below) and these variables pi can take 
values differing by many orders of magnitude (particu- 
larly in the solid phase) and must also satisfy the non- 
negativity constraint. This precludes the use of many 



standard minimization methods. The procedure we use 
here is the sam e . as . that originally employed in the hard 
sphere problemcacJ ! E3 with the important difference that 
the calculations involving the interacting term are per- 
formed in wavevector space. This is for two reasons: first, 
it turns out to be much more efficient, since the time used 
by Fourier transforming back and forth using efficient 
FFT routines turns out to be negligible compared with 
the time saved by not having to perform a double sum 
over the large computational lattice. Second, the direct 
correlation function is more conveniently computed in 
any case in terms of the azavevector variable. The proce- 
dure we use incorporates^ the nonnegativity constraints 
and is insensitive to the large range of the variables, but 
required 2 ^ a large number of iterations for convergence. 
The efficiency of the FFT method, however, still allows 
us to use the sizes required for the problem. 

The minimization procedure finds a local minimum of 
the free energy. The minimum located depends on the 
initial values chosen for the set of variables pi. The ap- 
propriate choices in each case are discussed in the next 
section. Generally speaking, nearly uniform initial con- 
ditions lead to liquid minima while ordered initial condi- 
tions with the proper symmetry lead to crystalline states. 

III. RESULTS 
A. General and one pin 

In this Section we present and discuss our numerical 
results. In principle, these could be given in terms of a 
minimal set of dimensionless parameters. However, it is 
more appropriate, in view of this paper's objectives as 
discussed in the Introduction, to present the results in 
terms of physical parameters with dimensions. This is 
the course we have taken. The values of the material pa- 
rameters that we use here have been therefore chosen as 
appropriate to BSCCO. These parameters are the pene- 
tration depth A and the interplane distance d, which to- 
gether with the temperature and fundamental constants 
determine T. We set d = 15.261, A(T = 0) = 1500A, 
and assume a two-fluid temperature dependence of A(T) 
with T C (H = 0) = 85K. For these parameter values, the 
dimensionless quantity T ~ 2650/ (T in Kelvins) at low 
temperatures where the T-dependence of A is negligible. 
We study temperatures and fields in the region where 
the melting transition of the vortex lattice is expected to 
occur. The strength of the pinning potential is described 
by the parameter a = Vq/T, as introduced above. We fix 
the range rg to = O.Iao, which corresponds to about 
55 A for B = 2kG. 

To test our procedures and to find out more about the 
parameter range to be studied and the system sizes re- 
quired, we begin by considering the simple case of an 
isolated pinning center in a vortex system in the liquid 
state. We place this pinning center in the center of the 
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FIG. 1. Numerical accuracy of density profile calculations. 
The normalized local density variable, p/po, in the presence 
of a single pinning center is plotted vs. r, the distance from 
the pinning center in units of ao. Two sets of data are plotted, 
for the same physical parameter values (B — 2kG, T = 20K). 
One set was computed with TV = 512, ao/h = 20 (crosses) and 
the other set (solid line) with TV = 1024, ao/h = 40 (see text). 
A third set with N = 2048, a /h = 80 would be completely 
obscured by the solid curve if plotted. 

computational lattice. Since periodic boundary condi- 
tions are used, this amounts to considering a periodic 
array of pinning centers with spacing equal to the size L 
of the computational cell. As discussed below, the values 
of L used in these studies are sufficiently large, so that 
the behavior near a pinning center is not affected by the 
presence of its periodic images. We then choose the initial 
configuration of the variables pi as one vortex located at 
the pinning center and uniform density everywhere else, 
with the average density consistent with po. We then 
perform the minimization of the free energy as described 
above. The main issues here are the determination of the 
appropriate values of the pinning strength, and, from a 
technical standpoint, finding the system sizes, and the 
values of h/ao required. A smaller value of h/ao implies 
higher spatial resolution in describing the variation of the 
local density, but at constant N this amounts to a reduc- 
tion in the size of the system being studied. One must 
therefore strike a balance. 

We have performed this procedure for fields B — 2kG 
and 3kG and at several temperatures in the range of in- 
terest (15 — 22K) at those fields. We have considered 
values of N ranging from 128 to 2048 with ao/h from a 
maximum of 80 down to values of order unity. Repre- 
sentative results of this rather extensive study are shown 
m Fig. |. In this Fi gure we plot the density variable p 
(normalized by po) &s a function of the dimensionless dis- 
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FIG. 2. Short distance (within the range of the pinning 
potential) results for the normalized density profile (dots). 
A fit to an exponential form in the pinning potential v p (r) 
(dotted line, see text) is also shown. 

tance r from the pin, measured in units of ao- The data 
in the Figure are all taken at a temperature T — 20. OK 
and a field B = 2kG, with the parameter a set at 0.06. 
Results are shown for two cases: N = 512 with ao/h = 20 
and N = 1024, a /h = 40. The scaling of a /h with N 
ensures that we are comparing systems containing the 
same number of vortices. The two results are very close 
to each other. We have also obtained results for N = 
2048, ao/h = 80, which are completely indistinguishable 
from those at N = 1024, so that if we had plotted them 
they would not be visible: the two plots at N = 1024 and 
N = 2048 would be on top of each other. This and simi- 
lar data obtained at other temperatures and fields tell us 
the range of values of N and ao/h needed to obtain high 
quality data. The results subsequently presented in this 
and the next subsection are obtained at N = 1024 and 
ao/h = 40 

The high peak at small r in Fig. [l] represents the large 
vortex density at the pinned site. This density then de- 
cays away in an oscillatory manner, as shown in the Fig- 
ure, towards its long range U*nit, which is unity for our 
normalization. As expectedlS, the behavior of p(r)/po 
outside the range of the pinning potential is very sim- 
ilar to that of the radial distribution functioned of the 
unpinned vortex liquid at the temperature and magnetic 
induction being considered. Thus, the medium and long 
range behavior of the density is reasonable. The behavior 
of p(r) at very short distances, inside the pinning range, 
is shown in Fig.0. One can see in this Figure how the 
results are well fit, as expected, by an exponential form 
e~ Vp ( r >, where the pinning potential in units of fc^T is 
given by Eqn. (p]q). This confirms that our computa- 
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FIG. 3. Integrated density within the range of the pinning 
potential (that is, the average number of vortices pinned at a 
center), as a function of pinning strength given by the param- 
eter a, defined in the text. For values of the pinning strength 
just before the kink in the curve, nearly one vortex is pinned 
at a center. The temperature is 20K and the field 2kG. 
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FIG. 4. Density profiles in the presence of two pinning cen- 
ters. The two centers are placed symmetrically around the 
origin (see text). In the main plot, the density profiles are 
shown as a function of the distance r from the origin, that 
is, from the midpoint between the two pinning sites. The in- 
set shows (solid line) the free energy as a function of the pin 
separation d. The horizontal dashed line marks the infinite 
distance limit. The curves shown in the main plot are for 
the cases where the distance d between centers equals 1.95ao 
(solid curve) and 2.60ao (dashed curve). The first corresponds 
to a minimum of @F (inset), the second to a nearby maximum. 
The results shown are for B = 2kG, T = 20K. 



tional mesh is sufficiently fine even at these very small 
ranges. 

Next, we must find the appropriate values of the pin- 
ning strength parameter a, as introduced above. We 
wish to consider here, of course, the case of strong pin- 
ning, but nevertheless a should remain small enough so 
that the total amount of flux pinned at each site remains 
on the average below one superconducting flux quantum. 
To choose the appropriate value, we studied the average 
number of vortices pinned as a function of a. Sample 
results are shown in Fig. ^, where the number of vortices 
pinned at the site (obtained by integrating the density 
over the pinning range) is displayed as a function of a. 
The data shown in this Figure were taken at B = 2kG 
and T = 18. OK; data at other relevant fields and temper- 
atures are very similar. We see that the number of pinned 
vortices rises very rapidly with a until a very marked kink 
occurs, at about the value when one vortex is pinned on 
the average. It is clear that one should use a value of a 
just below the kink in the curve, and in this study we 
have used the values a = 0.05, and a — 0.06 (the smaller- 
value was used in calculations at lower temperatures). 
For T in the range of interest here (15 — 22K), these val- 
ues of a correspond to Vq ~ 7 — 9. This is consistent 
with the results of the two-dimensional study of Ref. [l8| 
where it was found that the average number of vortices 
trapped at a pinning center decreases sharply below one 



as the dimensionless pinning strength Vq becomes lower 
than about 8. 



B. Two pinning centers 

Having in the previous subsection determined the 
properties of the density profile when the pins are, in ef- 
fect, very far apart, we consider now the case where there 
are two pinning centers separated by a smaller distance 
d. With a set so that each center pins nearly one vortex, 
we expect to have, when the two centers are not too far 
apart, interactive effects as a function of d, as the den- 
sity oscillations emanating from each pinning center (see 
Fig. must be distorted to match each other. We per- 
form this calculation by placing the two pinning centers 
symmetrically around the center of the computational 
triangular lattice, on the longer diagonal. For the initial 
conditions, we place one vortex on each of the two pinning 
sites, and a uniform density on the remaining computa- 
tional sites, consistent with the average density being po- 
We have performed this study at fields of 2kG and 3kG 
and at several temperatures. Results at B — 2kG and T 
— 20K, with a = 0.06 are shown in Fig. [|. These results 
were again obtained at N — 1024, do /h = 40. In the main 
plot we give two examples of the density profile as a func- 
tion of distance. Only one half of the density distribution 
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is plotted, with the origin corresponding to the center of 
the lattice and the horizontal axis representing distances 
along the diagonal. The distribution is then symmetric 
about the origin, and the distance from the center of the 
first peak to the plot's origin is half the interpin distance. 
The solid curve corresponds to the case where d (in units 
of ao) is 1.95 and the dashed curve to the case where it 
equals 2.60. One can see at larger values of r, away from 
both pins, an oscillatory decay similar to that in Fig. [I]. 
The behavior in the interpin region near the plot origin 
is more complicated. At shorter interpin distances (as in 
the solid curve), the density is very small between the 
two centers, but when that distance is increased to over 
twice ao (see the dashed curve), it becomes possible to 
have a peak between the two pinning centers, and as d 
is further increased, additional intersite peaks appear as 
well. 

When the two centers are close enough to interact, the 
free energy will obviously depend on whether the oscil- 
lations in the density profiles corresponding to the two 
centers "lock" or not. This impliesll3 that the free energy 
of the system should be an oscillatory function of inter- 
pin distance. To verify this, we have evaluated the free 
energy as a function of interpin distance d. Results are 
shown in the corner inset of Fig. ||. The dashed horizon- 
tal line is the result for the case where d is very large, 
that is, twice the value for a single pin (recall that the 
zero of free en ergy is taken to be the at the uniform liquid 
state, see Eq.($%). This value is 0F = -4.680 for the 
case plotted. The solid curve shows the behavior of (3F 
as a function of interpin distance and it clearly displays 
the oscillatory behavior of this quantity. As found in 
Ref. [l8], the free energy has minima at interpin distances 
approximately corresponding to the position of the max- 
ima of the single pin density profile shown in Fig. [l| This 
reflects that it is easier, for those distances, to lock the 
oscillations corresponding to the two centers. The two 
pin distances corresponding to the two curves shown in 
the main plot of Fig. ^ were chosen so that one corre- 
sponds to a free energy minimum (solid line) while the 
other (dashed curve) corresponds to a nearby maximum. 
The higher value of the free energy in the latter case is 
due to greater difficulty in matching the two profiles in 
this case. This difficulty is reflected in the smaller height 
of the first peak of p(r) outside the range of the pinning 
potential and the appearance of a small peak of p(r) near 
r/ao — 0.6. The oscillatory behavior of the free energy as 
a function of d implies an oscillatory dependence of the 
magnetization of the vortex liquid on the applied mag- 
netic field when a periodic array of pinning centers is 
present. In particular, the reversible magnetization in 
the liquid state is expected to show minima near certain 
integral values of BjB^. 

The integrated density inside the range of a pinning 
center remains close to unity when the pin separation d 
is greater than 2ao- As the value of d is decreased be- 
low 2ao, the integrated density begins to decrease and 
becomes substantially lower than unity for d < l.bao- 
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FIG. 5. Computation of the melting temperature for the 
pure system at a field B — 2kG. The symbols, connected 
by a solid line, represent the computational results for the 
dimensionless free energy of the crystal, as explained in the 
text. The results are shown to be independent of N which 
determines (see text) the mesh size used in the computation. 
The temperature at which the solid line crosses the liquid free 
energy (zero by convention, dashed line) is the melting point. 



Thus, the simultaneous occupation of two pining centers 
by two vortices is likely if the centers are far apart (in 
units of ao), but unlikely only when the two pinning cen- 
ters are separated by distances less than about 1.5ao- 
This is consistent with decoration experimental which 
show that nearly all pinning sites are occupied by vor- 
tices when the number of pinning sites is smaller than 
the number of vortices (B > B^,). 



C. Melting in the pure system 

As a preliminary step in our study of the effect of an 
array of pinning centers on vortex lattice melting, we 
have carried out calculations of the melting transition of 
the pure system for B = 2kG and 3kG. This is in order 
to determine the "clean" limit of our subsequent results. 
In addition, our numerical solutions in the pure limit 
can be coranaied with those obtained from variational 
treatmentsESEj of the same RY free-energy functional in 
which the density distribution in the crystalline phase 
was expressed in terms of the Fourier components of the 
density at a few small reciprocal lattice vectors. 

The computational cell used in our pure limit calcula- 
tions is one triangular-lattice unit cell with lattice con- 
stant L. The spacing h of the computational grid is cho- 
sen to have the values L/N with iV = 16, 32, 64 and 
128. Crystalline minima of the free energy are obtained 
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FIG. 6. Determination of the equilibrium lattice parameter 
Lo at melting. The dimensionless free energy of the crystal 
at fixed B and T is plotted here vs. the lattice parameter L. 
The free energy has a minimum at L = La ~ 1.988ao. 

by running the minimization routine with initial states 
that have a sharp peak of the density at the center of 
the computational cell. At sufficiently low temperatures 
such minima are found for a range of values of L. Typi- 
cal results obtained for B = 2kG, L = 1.9884ao, and two 
values of N (N = 16 and N = 128) are shown in Fig. | 
where the dimensionless free energy f3F of one unit cell 
of the vortex crystal is plotted as a function of the tem- 
perature T. We find that the free energies of the crystal 
obtained for all the values of N listed above are essen- 
tially the same, as exemplified by the data shown for two 
values of N, indicating that the effects of discretization 
are minimal provided that h < L/16 ~ 0.125ao. 

The equilibrium value Lq of the lattice parameter is 
determined by finding the value of L that minimizes the 
free energy at a given B and T. The dependence of /3F 
on the value of L for B = 2kG, T = 18. 5K is shown in 
Fig. [| The value of Lq is found to be close to 1.988ao, 

which is slightly higher than the spacing \J 2ir/\^3 ao 
of a triangular lattice of density po. This reflects the 
well-known resuhxl that the density of a vortex lattice in- 
creases slightly at melting. The transition temperature 
is determined from the zero-crossing of the free energy of 
the crystalline state, calculated for L — Lq, as a function 
of T, as illustrated in Fig. which shows the results of 
computations performed at L = Lq. For B = 2kG, the 
melting temperature is then T c = 18.45K. This value,-&f 
T c is slightly higher than that obtained variationallyc3. 
This is expected: the free energy of the crystal obtained 
in our unconstrained minimization should be lower than 
that obtained in calculations where the free energy is 
minimized with respect to a few parameters only. 



0.2 0.4 0.6 0.8 ' 

r/a 

FIG. 7. Radial dependence of the density distribution in 
the vortex lattice. The quantity plotted is the angular aver- 
age of the normalized local density p(r)/po. It is given as a 
function of the distance from the center of a crystalline unit 
cell. The symbols (crosses) are the computed results, and the 
solid line is a gaussian best fit, valid at small distances from 
the center of the cell. 

The discontinuity in the entropy at the crystallization 
transition is obtained from the numerically calculated 
slope of the (3F versus T curve at the transition temper- 
ature. The Fourier transform of the density distribution 
at the crystalline minimum obtained at the transition 
temperature yields the value of the jump in the crys- 
talline order parameter m, defined as the magnitude of 
the Fourier component of the density at the shortest re- 
ciprocal lattice vector of the triangular lattice. At B = 
2kG the entropy change As per vortex is 0.29/cs, and 
the jump in the order parameter m is Am = 0.52. Very 
similar results are obtained for B = 3kG: T c = 15.10K, 
As = 0.28fc B , Am = 0.52, L Q = 1.985a . Theseuaa; in 
close agreement with the results of earlier studiearcS. 

Our numerical method provides detailed and accurate 
information about the spatial distribution of the time av- 
eraged density in a unit cell of the vortex lattice. We find 
that the density function near the center of the unit cell 
is, to a good approximation, invariant under a rotation 
about an axis perpendicular to the layers and passing 
through the center of the unit cell. Fig. |?] shows a plot 
of the local density p(r) (angularly averaged and normal- 
ized by the average density po) at the crystalline mini- 
mum obtained for B = 2kG, T = 18. 5K, L = 1.988o , as 
a function of the distance r from the center of the unit 
cell. As shown in the figure, a gaussian fit to the data 
for small r provides a good account of the dependence 
of the density on r except at larger distances. The value 
of the Lindemann parameter C at melting may be ob- 
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tained approximately from the width of the gaussian fit, 
or more accurately from a numerical evaluation of the 
root-mean-square displacement 



/ r 2 p(r)dr 
J p{v)dv 



(3.1) 



where the integrals are over a lattice unit cell and r is 
the radius vector measured from the center of the unit 
cell. We find that the value of C at melting is 0.26 for B 
= 2kG and 0.25 for B = 3 Ml These values, similar to 
those found in earlier workBEj, are substantially larger 
than the typical values of £ in simple three-dimensional 
solids near melting. The large value of C implies that 
the peak of the density at the center of the unit cell is 
not very sharp (see Fig. [?]). This helps explain why rela- 
tively coarse values of the mesh size h (e.g. h = L /16) 
are adequate for obtaining an accurate description of the 
density distribution in the crystalline state. 



T C (K) 
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D. Periodic array of columnar pins 

Having obtained these clean limit results in the pre- 
vious subsection, we proceed now with our study of the 
effects of a commensurate, periodic array of pins on the 
vortex lattice melting transition. We consider a triangu- 
lar lattice of pins with spacing equal to ILq where I is 
an integer, and Lq, as defined above, is the equilibrium 
value of the spacing of the pure vortex lattice at its melt- 
ing point for the value of B being considered. Thus the 
pin concentration is c = I /I 2 . The computational cell 
used is one unit cell of the pin lattice (which contains I 2 
unit cells of the vortex lattice) with periodic boundary 
conditions and one pin located at the center of one of 
the vortex lattice unit cells. The value of h was fixed at 
Lo/64 in the calculations for B = 2kG. We also carried 
out a few calculations for B = 2kG using h = Lq/\Q. The 
results obtained for this larger value of h were found to 
be essentially the same as those obtained for h — Lo/64. 
We therefore used h = Lq/16 in our calculations for B 
— 3kG. We set the pin strength parameter, as mentioned 
above, to a = 0.06 in the calculations for B = 2kG. A 
slightly smaller value, a = 0.05, was used in the calcu- 
lations for B = 3kG because the temperature range of 
interest is lower in this case and the strength required to 
pin slightly less than one vortex is somewhat smaller. 

The crystalline and liquid minima of the free energy 
were located from "heating" and "cooling" runs, respec- 
tively. In a heating run, the crystalline minimum was 
first located by minimizing the free energy starting with 
an initial state in which the density distribution in each of 
the I 2 vortex lattice unit cells contained in the computa- 
tional cell was that in one unit cell of the vortex lattice of 
the pure system at its melting point for the same value of 
B. The crystalline minimum so obtained was "followed" 
to higher temperatures by increasing the temperature in 
small steps and running the minimization program at 



FIG. 8. Phase diagram (transition temperature T c as a 
function of pin concentration c) for B =2 kG (solid symbols) 
and B — 3kG (open symbols). The circles denote first order 
transitions and the squares mark crossovers. The dotted lines 
are polynomial fits, included to guide the eye. The arrow 
marks the approximate position of the critical point at B = 
2kG. At B — 3kG the critical point is near T = 17.6K 

each temperature with the minimum obtained at the pre- 
vious step as the initial state. In a cooling run, a liquid 
minimum was first obtained at a relatively high temper- 
ature by minimizing the free energy with an initial state 
consisting of one vortex located at the pinning center and 
uniform density everywhere else, so that the average den- 
sity was pq. This minimum was then "followed" to lower 
temperatures as in the heating runs, but decreasing the 
temperature in small steps instead of increasing it. The 
liquid state is not homogeneous in the presence of pinning 
and its free energy is nonzero. 

For both values of B studied and relatively small val- 
ues of c (I = 10, 8, 7 and 6), we found a range of tem- 
peratures over which both crystalline and liquid minima 
are locally stable. The crossing of the free energies of 
these two minima defines a first-order transition between 
crystalline and liquid states. Results for the transition 
temperature T c as a function of c for B = 2kG and 3kG 
are shown in Fig. |8[ The presence of columnar pins is 
found to increase T c . This should be expected: columnar 
pins suppress the disordering effects of the lateral wan- 
dering of vortex lines, and a commensurate periodic ar- 
ray of such pins clearly promotes crystallization. In other 
words, an external potential having the same symmetry 
as the crystal favors the crystalline state. The results for 
B = 2kG and 3kG are quite similar, with T c for B = 3kG 
reduced by approximately 3.4K for all these values of c. 
The discontinuities in the entropy s and the order param- 
eter m decrease as c increases!!!!! because pinning-induced 
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FIG. 9. Density distribution for the coexisting liquid and 
crystal states at the melting transition for pin concentration 
c = 1/49 and field B — 2kG. The normalized density profiles 
are plotted along a line joining two pins, one near x = 0, and 
the other just beyond the right edge of the Figure. 



FIG. 10. Sixfold angular modulation of the density in 
the liquid state for B = 2kG, T = 19.6K, c = 1/49. The 
normalized density averaged over the first peak (see text) is 
plotted vs. the angle measured from a line joining two pinning 
sites. 



order in the liquid increases with c. 

Our results yield not only bulk quantities but also very 
detailed information on the density distribution of the 
vortices. This quantity is experimentally accessible in 
scanning tunneling microscopy (STM) and scanning Hall 
probe measurements. In Fig. ^, we show the variation 
of the local density p along a line joining two neighbor- 
ing pinning centers for the crystalline and liquid minima 
near the transition temperature for c = 1/49 and B = 
2kG. The density profile in the liquid minimum can be 
viewed as a superposition of licjuid-like profiles near indi- 
vidual pins (compare with Fig. |l| noting the vertical axis 
scale), with small-amplitude, damped oscillations about 
p = Pq. In contrast, the density in the crystalline state 
exhibits higher, sharper and asymmetric peaks at the lat- 
tice points, with the density rising more sharply on the 
side closer to the pin, particularly at smaller distances 
from the pin site. Similar plots for c = 1/64 may be 
found in Rcf. 30 . These plots clearly bring out the obvi- 



ous differences between the structures of the coexisting 
crystal and inhomogeneous liquid phases. 

The density distributions at the liquid minima exhibit 
as expectecO a six-fold angular modulation. This is il- 
lustrated in Fig. |l0| where we have shown the average 
density at the first peak of the normalized local density 
near a pinning center (i.e. the density averaged over the 
region 1.75cto < r < 1.9ao where r is the distance from 
the pinning center) as a function of the angle measured 
from the line joining the pinning center to one of its near- 
est neighbors. The data shown are for the liquid mini- 
mum obtained for B = 2kG, T = 19. 6K and c = 1/49. A 



six-fold angular modulation of the density is clearly seen 
in the figure. The minima of the density occur on the 
lines that join the pinning center to its nearest neigh- 
bors. This is different from the behavior found in the 
crystalline minima where density maxima occur on the 
lines joining neighboring pinning sites (this can be seen 
at a different value of c from inspection of Fig. |l2| below). 

The full two-dimensional density distributions at the 
liquid and crystalline minima obtained near the transi- 
tion temperature for B = 2kG and c = 1/36 are shown 
as gray scale plots in Figs, [ll] and [lj, respectively. The 
plot for the liquid minimum exhibits the usual correla- 
tion "hole" around the pinning site at the center, and 
concentric "rings" of alternating high and low densities 
with six-fold angular modulation. The angular modula- 
tion at the first ring is less pronounced (and less obvi- 
ous in a gray scale plot) than that depicted in Fig. [l(] 
for c = 1/49. The plot for the crystalline minimum il- 
lustrates how the detailed structure of the periodically 
arranged crystalline density peaks changes with the dis- 
tance from the pinning site at the center. The ability 
of our numerical method to provide detailed information 
about the density distribution in highly inhomogeneous 
states is clearly illustrated in these figures as well as in 
Figs. I and 0. 

The degree of order in the liquid state increases with 
the pin concentration c, thereby decreasing the differ- 
ence between the crystalline and liquid minima. This 
has drastic consequences for the crystallization transi- 
tion. The behavior we find for c > 1/36 (I < 6) for both 
values of B is significantly different from that described 
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FIG. 11. Gray scale plot, as indicated, of the normalized 
density field p(r)/po at the liquid minimum for c = 1/36. A 
pinning center is located at the center of the shown unit cell 
of the pin lattice. 

above. For c = 1/25 and c = 1/16, the apparent min- 
ima obtained in heating and cooling runs have almost 
the same free energy, but somewhat different values of 
the order parameter m. We have shown in Fig. [l3] plots 
of m versus T obtained from heating and cooling runs for 
B = 2kG and c = 1/25 (circles and squares). The small 
difference in the heating and cooling values of m peaks 
at a temperature T = T x ~ 21.2K. In Fig. [l4] we show 
the two corresponding density profiles obtained at T = 
21. 2K. These plots are analogous to those in Fig. || for 
c = 1/49. In sharp contrast to that case, (and also to the 
c = 1/36 case shown in Figs, [□] and |lj) the two profiles 
are now very similar, with the one obtained in the heat- 
ing run exhibiting only a slightly higher degree of order, 
consistent with the higher value of m. This leads to the 
suspicion that the free energy at c = 1/25 may have only 
one very "flat" minimum near T — T x under the con- 
ditions studied. When attempting to find a minimum, 
our numerical routine stops when the free energy gra- 
dient becomes smaller than a certain small convergence 
parameter. When a minimum is very flat, it may stop 
at slightly different configurations when approaching it 
from different directions in configuration space. 

If this situation occurs, then the density configuration 
at the true minimum of the free energy should be bet- 
ter approximated by a linear combination of the density 
configurations found in heating and cooling runs. We 
therefore evaluated the free energy for a set of configura- 
tions {pi(x)} defined by 

Pi {x) = xpV+{l-x)pf\ (3.2) 

where and {p^} are the density configurations, at 



FIG. 12. Gray scale plot of the normalized density field 
at the crystal minimum coexisting with the liquid minimum 
shown in the previous Figure. 

the same temperature, at the apparent minima obtained 
in heating and cooling runs, respectively. The mixing 
parameter x is in the range < x < 1. If one then plots 
the free energy thus obtained either as a function of x 
or, equivalently, as a function of m{x) = im' 1 ' + (1 — 
x)m( 2 \ where rrS^ 1 and mP^ are the order parameters 
in the two configurations, one finds that it indeed has 
a minimum at x — xq ~ 0.5 at temperatures near the 
expected transition, for all higher concentrations, c > 
1/25. An illustration, for the case B = 2kG, T = 21. 2K, 
and c = 1/25 discussed above, is provided by Fig. |l5|. If 
one attempts such a procedure at lower concentrations, 
on the other hand, the resulting plot turns out to have a 
maximum, rather than a minimum, near x = 0.5. 

This analysis shows that the suspicions mentioned 
above were correct and that for c > 1/25, only one mini- 
mum of the free energy exists at each temperature. The 
value of the free energy at this minimum is lower than 
those found in the heating and cooling runs. Thus, there 
is no first-order transition at c = 1/25 or higher. The 
line of first-order transitions found for smaller values of c 
ends at a critical point which lies between c = 1 /36 and 
c = 1/25 at both values of the field considered. 

At c > 1/36, above the critical point, a crossover 
rather than a sharp transition characterizes the change 
from liquid-like to solid-like behavior. The crossover tem- 
perature can be conveniently defined from the numeri- 
cally calculated temperature derivative of the "equilib- 
rium" value, m(xo), of the order parameter. Plots of 
both m (xg) and its temperature derivative are shown 
in Fig. [l|. The temperature at which the derivative of 
the order parameter peaks is obviously very close to the 
temperature T x defined earlier as that at which the dif- 
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FIG. 13. Illustration of the different values of the order 
parameter m found in heating and cooling runs (squares and 
circles, respectively) at higher values of the pin concentration 
(c = 1/25, B — 2kG in this case). The triangles represent 
the equilibrium values of m found as explained in the text. 
The solid line is a polynomial fit to the equilibrium data. The 
dotted curve is the absolute value of the derivative of the solid 
line. 



ference between the order parameters obtained in heat- 
ing and cooling runs peaks. The crossover temperature 
can therefore be identified with T x . The sharpness of the 
crossover suggests that c = 1/25, T = T x ~ 21.2K is close 
to the critical point for B — 2kG, as indicated by the ar- 
row in Fig. |^. For c = 1/16 the crossover is smoother. 
Our results for B = 3kG are very similar to those at the 
lower field, with a similar value of the critical c but lower 
crossover temperatures (this is obvious from Fig. |^), with 
T x ~ 17.6K for c= 1/25. 

The crossover to the crystal state at c > 1/36 may be 
correlated with the onset of percolation of vortices which 
are "localized" according to a density-based criterion. 
Localization of vortices is, strictly speaking, a dynam- 
ical phenomenon, but some information about localiza- 
tion may be obtained from the distribution of the time- 
averaged local density. The local density near a point 
where a vortex is localized should be significantly higher 
than the average density po- Therefore, a density-based 
criterion for localization may be obtained by demanding 
that the density near a point where a vortex is localized 
exceed a suitably chosen cutoff value p c . Our results for 
the density modulation around an isolated pinning cen- 
ter suggest an appropriate choice for this cutoff. One can 
see, for example, in Fig. [IJ that in the temperature range 
of interest the local density in the neighborhood of an 
isolated pinning center (but outside the range of its pin- 
ning potential) does not exceed 3po if the system is in the 



FIG. 14. Density profiles, depicted as in Fig. for the 
heating and cooling runs shown in Fig. [l3|, at temperature T 
= 21.2K, which is very close to the crossover temperature T x 
defined in the text. 



liquid state. This suggests that values of the local den- 
sity p less than 3po correspond to mobile vortices. We, 
therefore, take p c = 3po- We divide the computational 
cell into I 2 vortex-lattice unit cells and associate a local- 
ized vortex with a unit cell if the local density exceeds p c 
at some point inside that cell. We then examine whether 
the unit cells that contain localized vortices according to 
this criterion percolate across the sample. 

All the vortex lattice unit cells in a crystalline mini- 
mum contain localized vortices, since the maximum value 
of pi at the lattice sites of a crystal always exceeds p c . 
In contrast, only the vortex lattice unit cells in which 
pinning centers are located and, in some cases, the near- 
est neighbors of such unit cells, contain localized vortices 
in the coexisting liquid minimum at the crystallization 
transition for c < 1/36. Thus, for c < 1/36, the crystal- 
lization transition trivially coincides with a percolation 
of unit cells containing localized vortices. For c > 1/25, 
in the crossover region, we have found that the unit cells 
containing localized vortices do not percolate if the tem- 
perature is higher than the crossover temperature T x de- 
fined above, but percolation occurs below T x . Typical 
results are shown in Figs. 16 and fL7[ In Fig. [lq, we have 



shown the locations of the units cells containing localized 
vortices in the minimum obtained in the heating run for 
B = 2kG, T = 21. 2K and c = 1/25. These unit cells do 
not percolate across the 5x5 computational cell, while 
the cells containing mobile vortices do. Since the degree 
of localization in the minimum obtained in the heating 
run is higher than that in the minimum obtained in the 
cooling run, no percolation would be obtained at this 
temperature if the cooling-run minimum or the cquilib- 
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FIG. 15. The "mixed" free energy plotted as a function of 
m(x), as explained in the text, at c = 1/25, B = 2kG, T = 
21. 2K. The triangles represent the results of the computation. 
The slight irregularity of the data points reflects numerical 
uncertainties. The dashed line is a fit to the Landau expansion 
of Eq. (p|). 



rium configuration pi(xo) were used for finding the unit 
cells containing localized vortices. On the other hand, as 
shown in Fig. |lj], the unit cells containing localized vor- 
tices do percolate across the sample in the "equilibrium" 
configuration obtained at the slightly lower T — 21. IK. 
At this temperature, the heating run shows percolation, 
but the cooling run does not. Thus, percolation occurs 
at a temperature very close to the crossover temperature 
T x ~ 21. 2K. Very similar results were obtained for B 
= 2kG, c = 1/16, and B = 3kG, c = 1/25, indicating 
that this is a general condition. This result is physically 
reasonable: a system in which localized vortices perco- 
late (and consequently, the mobile ones do not percolate) 
should behave as a "solid" at long length scales. 

It is easy to see that the occurrence of a critical point 
in the phase diagrams of Fig. || does not contradict any 
fundamental principles. In the presence of commensu- 
rate periodic pinning, the liquid and the crystal have the 
same symmetry. Since the degree of order in the liquid 
increases with c, it is possible for the liquid and the crys- 
tal to become indistinguishable beyond a critical value of 
c. One can then go from one phase to the other without 
crossing a sharp phase boundary. 

The basic features of the phase diagrams may be un- 
derstood from a simple. Landau theory. From well-known 
symmetry argumentsEJ one can write down a Landau ex- 
pansion for F: 



BF = — a 7 m 2 a^m 3 H — 04m 4 

2 3 4 



gm, 



(3.3) 



where m is our order parameter, the constants 03 and 



FIG. 16. Location of cells containing localized vortices un- 
der the indicated conditions. The circles denote the positions 
of the 25 vortex lattice unit cells contained in a unit cell of the 
pin lattice. A pinning center is located in the unit cell at the 
bottom left corner. A star in a circle denotes that the unit cell 
contains a "localized" vortex according to the criterion given 
in the text. The temperature here is slightly higher than the 
crossover temperature T x for B = 2kG, c = 1/25. The cells 
containing localized vortices do not percolate across the 5x5 
sample. 



(I4 are positive, and ei2 decreases with T. The "order- 
ing field" g is proportional to the pin concentration c. 
A simple analysis shows that this free energy leads to a 
first-order transition for g < g c = a\jTJa\ and a criti- 
cal point at g = g c , a-i = a-2 C = a 2 /3ci4. The transition 
temperature increases with the ordering field g in agree- 
ment with the arguments previously discussed. The la- 
tent heat and the order parameter discontinuity Am van- 
ish as (g c — g) 1 / 2 as g approaches g c from below. It was 
shown in Ref. ^ojthat our data for As and Am are indeed 
well-described by the form oc (c c — c) 1 / 2 with c c close to 
1/25. 

More quantitatively, it is possible to fit our j3F vs. 
m data to the form Eqn. ([3/J). The best fit for B = 
2kG, T = 21. 2K, c = 1/25, is shown in Fig. [H]. The 
fitting parameter values are 0,2 = 70.71, 03 = 234.87, 
a4 = 261.71, and g = 7.12. Using the values of 03 and 04 
obtained from the fit, we get a^ c — 70.26, and g c = 7.01. 
These values are very close to, but slightly lower than 
the best-fit values of ai and g, indicating that the critical 
point for B — 2kG is, as we had already stated, very close 
to c = 1/25, T = 21. 2K. This explains the sharpness of 
the crossover at c = 1/25. The numerical results for B 
= 3kG can be analyzed in the same way. A fit for T = 
17. 6K, c = 1/25 yields then values of a^c and g c which 
are less than 1% lower than the best-fit values of 0,2 and 
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3=2kG, T=21.1K c=1/25 
* "localized" vortices 



o m o 



"equilibrium" 



FIG. 17. As in the previous Figure, but at a temperature 
slightly lower than the crossover temperature T x . The cells 
containing localized vortices now percolate across the sample. 



g, respectively. Therefore the Landau free energy gives 
a good semi-quantitative account of the critical behavior 
of our density functional computations for both values of 
B. This strengthens our conclusions about the existence 
and location of the critical point. 



IV. SUMMARY AND DISCUSSION 

We have used in this paper numerical minimization 
of a discretized free energy functional to study the ef- 
fects of columnar pinning on the structure and thermo- 
dynamics of a system of pancake vortices in the mixed 
phase of highly anisotropic layered superconductors. The 
most salient result of our study is the existence of a crit- 
ical point in the vortex-lattice melting phase diagram 
when a commensurate, periodic array of pinning centers 
is present. Our results show that the line of the melting 
transition in the T—c plane, which is of course first-order 
for small values of the concentration of pinning centers, 
terminates at a critical point as the pin concentration is 
increased. Beyond this critical point, the transition is 
replaced by a crossover, with smooth behavior of the or- 
der parameter and other thermodynamic quantities. To 
our knowledge, this is the first quantitative theoretical 
prediction of a continuous melting transition in a three- 
dimensional system. This critical point should be experi- 
mentally accessible: the pin lattice spacing for B = 2kG, 
c = 1/25 should be ~ 0.55 /im, close to the spacing of the 
radiation-induced, pin array of Ref. [l| The same group 
recently showedEj that columnar pins can be created in 
BSCCO in a highly controlled manner. We, therefore, 
expect that the fabrication of bulk HTSC samples with a 



periodic array of columnar pins with appropriate spacing 
is technically feasible. Fabrication of such samples and 
experiments to verify our theoretical predictions would 
be most welcome. 

We have shown that most of the features of our phase 
diagram can be understood from a simple Landau the- 
ory. The critical point found in our study is analogous to 
the liquid-gas critical point in mean-field theory. Fluc- 
tuations are expected to change this correspondence be- 
cause the symmetry of our order parameter is different 
from that for the liquid-gas transition. Theoretical stud- 
ies of the universality class of this critical point would be 
interesting. However, the location of the first order melt- 
ing line and the existence and experimental accessibility 
of the critical point should be quantitatively described by 
our work. We have also shown that the smooth crossover 
from liquid-like to solid-like behavior beyond the critical 
point might be interpreted as a percolation threshold for 
localized vortices. 

The one-pin results reported here provide useful infor- 
mation about the dependence of the average number of 
vortices trapped at a pinning center on the temperature 
and the strength of the pinning potential, while our re- 
sults for the interaction between two neighboring pins il- 
lustrate the occurrence of interesting effects in the liquid 
state arising from the commensurability of the separa- 
tion of the pins with the average inter-vortex separation. 
Finally, our method yields very detailed results for the 
density distribution in the system, which is accessible 
through STM and scanning Hall probe measurements. 

As noted in Section [0], the symmetry of the three- 
dimensional system considered here makes the calcula- 
tions effectively two-dimensional, with the function C 
playing the role of the two-dimensional direct pair cor- 
relation function. The direct pak-cerrelation function of 
a two-dimensional vortex Iiquidll3t3 is quite similar to 
the function C used in the present study. We, therefore, 
expect that most of the results obtained here should ap- 
ply, at least qualitatively, to thin-film superconductors 
in the presence of strong pinning centers. As noted in 
Section III , some of our one- and two-pin results are in- 
deed in good quantitative agreement with those obtained 
in Ref. |l8| for a two-dimensional system of vortices with 
strong pinning. This leads us to expect that the phase di- 
agram obtained here for the vortex lattice melting transi- 
tion in the presence of a periodic array of pins should ap- 
ply, with no more than fairly minor quantitative changes, 
to the melting transition of a two-dimensional vortex lat- 
tice in thin-film superconductors with commensurate pe- 
riodic pinning. Since periodic arrajzs, pf, strong pinning 
centers have already been fabricatedt3~tll in thin-film su- 
perconductors, our predictions can be readily tested in 
experiments. One should, however, keep in mind that the 
predictions of our mean-field-like density functional cal- 
culation are less reliable in two dimensions where the ef- 
fects of fluctuations are stronger. The melting transition 
in a pure twQ-dimcnsional system without pinning can be 
continuously, whereas the mean-field prediction of first- 
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order melting is always realized in pure three-dimensional 
systems. Our main result about the existence of a criti- 
cal point in the phase diagram should apply to thin-film 
superconductors with commensurate periodic pinning if 
the system exhibits a first-order melting transition in the 
absence of pinning. 

The melting transition of the lattice of interstitial vor- 
tices in the presence of a commensurate, periodic array 
of pinning centers provides a physical example of melting 
in the presence of an external periodic potential. Simi- 
lar melting transitions are of interest in other systems 
such as atoms adsorbed on crystalline substratesEH, and 
colloidal particles in interfering laser fields^ and arrays 
of optical trapsSl. Our method and results would be of 
relevance to these systems also. 
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